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Abstract 


This paper presents a method for the detection of traveling waves of activity in neural recordings 
from multielectrode arrays. The method converts local field potential measurements into the phase 
domain and fits a series of linear models to find planar traveling waves of activity. Here I present 
the new approach in the context of the previous work it extends, apply the approach to data from 
neural recordings from a single animal, and verify the success of the method on simulated data. 


1 Background and Overview 


Current technology allows for neural recording 
on many scales, each with trade offs in terms 
of spatial resolution, temporal resolution, and 
recording target. A single extracellular electrode 
has high temporal resolution and relatively high 
spatial precision. Electrode arrays, by recording 
from spatially structured collections of electrodes 
simultaneously, supplement this with some infor¬ 
mation on spatio-temporal patterns. Electrode 
arrays may therefore be a useful recording tech¬ 
nique for exploring local network dynamics in 
vivo. 


The signal recorded by a single electrode 
contains high frequency electrical signals from 
nearby spiking neurons and a low frequency sig¬ 
nal referred to as the local field potential (LFP). 
The LFP is believed to be a measure of the 
summed synaptic activity in a radius around 
an electrode (iKatzner et al.l 120091). Previous 


work in the Hatsopoulos lab at the University 


of Chicago (Rubino et ah, 2006) suggests that 


LFP recorded over an array of electrodes can re¬ 
veal spatio-temporal patterns in activity. Specif¬ 
ically, it can reveal traveling waves of oscillation. 
Traveling waves are widely discussed in the 


neuroscience literature (Delaney et al. 


1994 


Huang et ah, 2004[ Manjarrez et al.l |20071 |Wu 

et al., 

20081 ^1-' 

20071 Reimer et ah, 2010). 


Propagating waves of oscillation have been re¬ 
ported in the visual, auditory, olfactory and so¬ 
matosensory systems as well as in motor cortex, 
hippocampus, and spontaneous cortical activity. 


Theoretical work reported in Ermentrout and 


Kleinfeld (2001) demonstrates that the traveling 


wave patterns reported in experimental studies 
are likely to arise in realistic neural networks. 
Even so, the possible computational purpose, if 
any, of traveling waves remains unclear. Deter¬ 
mining their significance, or even relegating them 
to the realm of the epiphenomenal, will require 
advances in the methodology used to detect and 
describe them. 

Currently, there is no generally established 
methodology for detecting and describing trav¬ 
eling waves. Many studies rely on visual in¬ 
spection of Voltage Sensitive Dye (VSD) im¬ 
ages. Several recent papers using LFP electrode 
recordings make a start on developing formal 
wave detection techniques. [Lubenov and Sia^ 


pas (2009) used regression methods on a collec¬ 


tion of electrodes, and Rubino et al. (2006) used 


phase derivative calculations on electrode array 
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data. This project synthesizes ideas from these 
two papers (applied to the data from the later) 
to develop methodology that allows for statisti¬ 
cally rigorous single trial analysis of LFP array 
data. 

This research analyzes data from the Hat- 
sopoulos Lab at the University of Chicago. I 
first replicated the results reported in |Rnbino 


et al. (2006) regarding the detection of travel¬ 


ing waves. This previous methodology was able 
to detect the presence of wavelike activity and 
capture its broad characteristics, but was unable 
to provide a single trial description of activity. 
I made several improvements to the methodol¬ 
ogy, allowing for single trial analysis. I also 
demonstrated that wavelike activity is persistent 
rather than occurring in short isolated time pe¬ 
riods as previously reported. Finally, I created 
artificial data to test both the original and im¬ 
proved methodology under varying noise condi¬ 
tions, demonstrating the efficacy of my improve¬ 
ments. 


The next two sections of this paper are intro¬ 
ductory. Section Instates the model for a planar 
wave of oscillation. Section |3] describes the data 
being modeled. The following section discusses 
the results of previous analysis, while section 
describes my analysis and results in detail. Fi¬ 
nally, section describes the conclusions from 
artificial data analysis. 


2 Traveling Waves 


A traveling wave is composed of a collection of 
oscillators whose relative phases depend on their 
spatial location. The simplest form of this phase 
relationship is a planar wave, where the phase of 
the oscillators is linear in space and time. With¬ 
out loss of generality, we can set the phase at 
time 0 and position 0 to be 0. In one dimension, 
the equation to describe a planar wave can then 
be written 


V{x^t) oc sin 


27r 



( 1 ) 


where / is the frequency of the oscillators and s 
is the speed of the wave. 

This linear phase relationship can be easily ex¬ 
tended to describe planar waves in two dimen¬ 
sions by making x a vector. With two spatial 
dimensions, more complex patters are possible 
such as bull’s eye or spiral waves. However, this 
research addresses only the simpler, linear case. 

Since the phase is linear, fitting this model re¬ 
quires estimating three parameters: the change 
in phase in the x, y and t dimensions. From 
these parameters we can calculate properties of 
the wave that are of scientific interest such as its 
speed and direction. Table shows this relation¬ 
ship. 


3 Data 


The data I analyze here were recorded by Nicho 
Hatsopoulos’ lab at the University of Chicago. 
They consist of 391 trials of a center-out reach¬ 
ing task performed by a monkey. On each trial, 
the monkey was given a target in one of eight di¬ 
rections (the “instruction”) and then had to wait 
for one second before initiating movement to the 
target (after the “go” cue). 

LFP voltages were recorded on a 96 channel 
Utah array implanted in motor cortex. The data 
were sampled at 1 kHz. Based on prior research 
(Rubino et ah, 2006) I chose the waiting period 
between the “instruction” and “go” cues to be 
the period of interest for further analysis. As 
shown in Figure this period is characterized 
by strong oscillatory activity in the Beta band 
between 15 and 20 Hz. 

Fourier analysis of the waiting period data re¬ 
vealed the expected peak in the Beta band. Ad¬ 
ditionally, the Fourier spectrum revealed signifi¬ 
cant high frequency noise at multiples of 60 Hz, 
especially in two channels. This noise is pre¬ 
sumably electrical line noise, and is removed by 
filtering during analysis. 

The 96 channels on the array are very highly 
correlated. Figure shows the voltages for all 
channels superimposed for one trial. For this 
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Parameter 

Interpretation 

Calculation 

Frequency 

The frequency of the oscillators composing the wave. The 
current model requires all oscillators to have the same fre¬ 
quency. 

/ OC (5t 

Direction 

Wave propagation direction. 

d oc atan 

Speed 

Spatial speed of contours of constant phase. Higher speeds 
indicate more synchronous activity. 

Sx -\-Sy 


Table 1: Calculating wave parameters from linear model parameters. Here 6x, 
5y, and 6t are the change in phase in the x, y, and t dimensions. Wave speed 
and direction are the primary parameters of scientific interest. 
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Fig. 1 : This spectrogram, computed using Morlet wavelets, shows high power 
between 15 and 20 Hz during the waiting period (between 2000 and 3000 ms in 
this figure). The spectrogram was computed individually for each trial and then 
averaged over all trails, aligning on the time at which the monkey received the 
direction cue. 
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figure, the data were wide-band filtered between 
10 and 45 Hz for clarity, but the correlation of the 
channels and strong oscillatory activity around 
20 Hz is quite clear. 


4 Summary of Previous Work 


Previous work reported in Rubino et al. (2006) 


found evidence of wavelike activity. The direc¬ 
tions of wave propagation were seen to follow a 
bimodal distribution with peaks approximately 
TT radians apart (nearly opposite directions). The 
wavelike activity was found to consist of short 
periods of waves surrounded by non-wavelike pe¬ 
riods. Here I give a very brief summary of these 
results and the analysis used to reach them. 
More details can be found in the Analysis sec¬ 
tion or in the original paper. 

As described in sectionfitting a planar wave 
model to data requires estimating three param¬ 
eters: 6x^ Sy^ and 6t (the change in phase in each 
of the three dimensions). The first step of the 
analysis is to extract the phase. This is done 
using the analytic signal, discussed in detail in 
section [5T| At each point in space and time, the 
local gradient of the phase is computed. These 
gradients are then averaged spatially to give es¬ 
timates 6x^ Sy^ and 6t at each time point. 

In addition to fitting the three parameters, it 
is necessary to evaluate the model fit and deter¬ 
mine whether activity is indeed wavelike. Previ¬ 
ous work used a measure of phase gradient align¬ 
ment. If the data are actually a planar wave, 
then the phase will be linear in space and time 
and the spatial phase gradients at all locations 
will point in the same direction. If the data are 
not wavelike, the spatial gradients would be un¬ 
likely to align. Figure 1^ demonstrates this con¬ 
trast. 

To measure phase alignment, Rubino et al. 
used a measure they called Phase Gradient Di¬ 
rectionality (PGD). 


PGD = 


\n '^1=1 






( 2 ) 



Fig. 3: Two time points demonstrating high 
and low phase gradient alignment. Eaeh 
hlaek line indieates the direetion of the spa¬ 
tial gradient on one of the ehannels of the 
array at the time point. The red line is the 
mean over the array. The length of the red 
line is the MRL, a measure of phase gra¬ 
dient alignment. At time 1 the gradients 
are well aligned (indieating wavelike aetiv- 
ity). At time 2, the gradients are not aligned 
(indieating non wavelike aetivity) 


where n ranges over electrodes. This is very sim¬ 
ilar to the standard circular variance measure 


known as Mean Resultant Length or MRL (Mar 


dia and Jupp, 1999). 




Wy^xUOyi) 


( 3 ) 


The only difference is that PGD weights individ¬ 
ual gradients according to their magnitude while 
MRL normalizes all gradients. The two measures 
are very similar for practical purposes. 

The alignment measure is computed for each 
time point. Since wavelike activity should have 
high alignment, Rubino et al. set a threshold of 
0.5 and classified times at which the PGD was 
above this value to be wavelike. The alignment 
measure for an example trial is shown in figure 

SI 

Once the model is fit to the data, it is possible 
to examine the wave characteristics such as prop¬ 
agation direction. Figure shows histograms 
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Fig. 2: Voltage from all channels superimposed for one trial The data were wide¬ 
band filtered between 10 and 45 Hz before being plotted. 



Fig. 4: Alignment measures computed for an example trial. Previous work 
identifies waves and time periods when the alignment measure is above 0.5. 
This leads to the identification of short periods of wavelike activity, such as 
those highlighted in the figure. The alignment measures appear quite noisy, 
making single trial interpretation questionable. 
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of wave propagation direction for all times and 
for times classified as waves. Especially during 
the wavelike periods, the propagation directions 
show a clear bimodal pattern. 

The conclusions from previous work were that 
the data show short periods of wavelike activ¬ 
ity during which the waves propagate in two 
primary directions. I performed a permutation 
test to determine if the PGD threshold of 0.5 
was appropriate for wave identification. By per¬ 
muting the spatial organization of the channels 
I destroyed any spatial wavelike pattern while 
maintaining the temporal structure of the indi¬ 
vidual channels and the overall level of correla¬ 
tion across the array. In this way I computed 
a null distribution for the alignment measures, 
shown in figure 

Comparing the true PGD distribution to this 
null distribution, I concluded that the thresh¬ 
old of 0.5 was inappropriate. A better threshold 
would be closer to 0.2. This indicates that the 
wavelike activity in the data is a more perva¬ 
sive and consistent pattern in the system than 
previous work suggested. The rest of this paper 
describes the analysis methods I used to extract 
interpretable single trial results showing nearly 
constant wavelike activity. 


5 Improved Analysis Methods 
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A planar traveling wave of oscillation is de¬ 
fined by linear dependence of phase on spatio- 
temporal position. For this reason, a natural ap¬ 
proach to detecting such waves is to work in the 
phase domain. This approach is taken in papers 
such as iRubino et al. (2006) and Lubenov and 


Siapas (2009). The first step of this modeling 


approach is to extract the phase at each point in 
space and time. This process is described in sec¬ 
tion [Q Once the data is transformed into the 
phase domain, the remaining analysis involves 
fitting a linear model and evaluating model fit. 


Fig. 5: Histograms of estimated wave prop¬ 
agation direetion. (A) shows estimated di- 
reetions over all times and trials. (B) shows 
estimated direetions for those times identi¬ 
fied as being wavelike. The bimodal pattern 
of propagation direetions is elear, espeeially 
in B. 
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Fig. 6: The results of a permutation proeedure to derive a null distribution for 
the alignment measure (PGD) used in previous work. The grey distribution is 
the null distribution, and the blaek line indieates the 99th pereentile of the null 
distribution. This indieates that an appropriate threshold for the PGD would 
be 0.2 rather than 0.5. 


5.1 Extracting Phase 

The first step of a phase-based analysis is to 
extract the instantaneous phase at each spatio- 
temporal point. This is typically accomplished 
using the Hilbert transform to calculate the an¬ 
alytic signal, as described below. 

The Fourier transform of a real-valued func¬ 
tion is Hermitian. Because of this, the infor¬ 
mation about negative frequencies is redundant. 
The analytic signal discards these negative fre¬ 
quency components at the cost of converting the 
associated signal from real-valued to complex¬ 
valued. If x{t) is a real valued function with 
Fourier spectrum X(/), the corresponding ana¬ 
lytic signal is defined to be 

Xa{t) = x{t) + ix{t) (4) 

where x{t) is the Hilbert transform of x{t). 
The Hilbert transform is the convolution of x{t) 
with We can therefore see that equation 0 
becomes 


Xa{t) = x{t) + i (^x{t) * Pj (5) 

= x{t) * (6) 

Taking the Fourier transform of this, and 
defining h{f) to be the Heaviside step function, 
we see that the negative frequencies are 0, as ex¬ 
pected. 

Xa{f)^X{f)-2h{f) (7) 

2Xif) f>0 
X{f) / = 0 (8) 

0 /<0 

Returning to the time domain analytic signal, 
we can write equation @ in polar coordinates as 

Xa{t) = (9) 

where A{t) is called the amplitude envelope 
and (j){t) is the instantaneous phase, wrapped 


7 












Signal and A(t) 



Time 


Instantaneous Phase 



“I-1-1-1-1-1— 

0 200 400 600 800 1000 


Time 

Fig. 7: The instantaneous amplitude and 
phase ealeulated via the analytie signal for 
an amplitude modulated pure sine wave. The 
upper graph shows the original signal in hlaek 
with the amplitude envelope A{t) plotted in 
red. The lower graph shows the instanta¬ 
neous wrapped phase (l){t). 

between —tt and tt. Figure shows the instan¬ 
taneous amplitude and phase described by the 
analytic signal for an amplitude modulated pure 
sine wave. 

One concern in using the analytic signal to ex¬ 
tract phase is that is only interpretable as in¬ 
stantaneous phase for signals with narrow band 
frequency content. The data in this research is 
not a pure sine wave, so this characteristic of the 
analytic signal is a potential problem. In addi¬ 
tion, the model for a planar wave has only one 
frequency term, requiring a single frequency for 
all oscillators. However, since the data do show 
a well-defined peak frequency, it is reasonable 
to filter the data and look for wavelike activity 
around that peak. Previous work used a wide¬ 
band filter from 10-45Hz. My work shows that 
the previous filter was too wide. Interpretable 
results require a much narrower band filter; the 


results reported here used a filter with a 1 Hz 
passband (plus rolloff). 

5.2 Unwrapping Phase 

As shown in the previous section, the phase as 
extracted by the analytic signal is wrapped be¬ 
tween —TT and TT. Fitting a model to the data re¬ 
quires removing the 2tt jumps in the phase. This 
is a well-studied problem known as “phase un¬ 
wrapping” . With data satisfying strong assump¬ 
tions about the sampling rate and noise levels, 
phase unwrapping is a trivial problem. However, 
with realistic data, the problem is unidentifiable. 
The multitudes of unwrapping algorithms that 
have been proposed aim to find a good solution 
in the realistic case. 

To unambiguously unwrap the phase, we must 
make the assumption that the data were sampled 
according to the Nyquist-Shannon sampling the¬ 
orem in all dimensions; that is, that the data 
were sampled at a minimum rate of twice the 
highest frequency present in the data. If this as¬ 
sumption is satisfied, the true phase difference 
between any adjacent points will be less that tt. 
Therefore, unwrapping can be done by travers¬ 
ing the data with some spanning tree, adding or 
subtracting units of 2tt to a point if its value dif¬ 
fers from that of the previous point by more than 

TT. 

When local noise is introduced into the data, it 
may cause violations of the assumption. The ba¬ 
sic unwrapping algorithm may erroneously add 
or subtract 2tt or fail to do so when it should. 
Since the unwrapping is done along some tree, 
any points that are unwrapped subsequent to 
the error along the same branch of the tree will 
have the same error. This error propagation can 
result in wildly inaccurate unwrapping results. 
The general approach of any improved unwrap¬ 
ping algorithm is to chose the structure of the 
unwrapping tree such that noisy areas are un¬ 
wrapped last and error is not propagated. 

Previous work used very simple phase unwrap¬ 
ping, repeatedly unwrapping small portions of 
the data in one dimension at a time. As de- 
















scribed in section]^ Rubino et al. (2006) esti¬ 
mated the model parameters by computing the 
local phase gradient at each point and then av¬ 
eraging over the array. Computing each compo¬ 
nent of the gradient only requires using several 
data points in a line in a single dimension and 
so it is only strictly necessary to unwrap those 
points. Rubino et al. unwrap the necessary 
data points, compute the gradient component, 
discard the unwrapping results and then repeat 
the process for each dimension at each point. 
This unwrapping technique restricts dramatic er¬ 
ror propagation in the unwrapping by only un¬ 
wrapping small portions of data at a time. How¬ 
ever, it does nothing to limit propagation of er¬ 
ror within the local regions and restricts possible 
model fitting techniques to those that do not re¬ 
quire simultaneous consideration of a large set of 
data points (precluding, for instance, regression 
approaches). 


It can be demonstrated that considering all 
available dimensions during unwrapping can 
yield better results. This is because unwrap¬ 
ping requires choosing a path or tree to traverse 
the data. In one dimension, the path is prede¬ 
termined. In higher dimensions, there are in¬ 
creasing numbers of paths to choose from, in¬ 
creasing the likelihood that one can be found 
that avoids traveling through noisy regions of 
data. There are dozens of unwrapping algo¬ 
rithms that have been proposed for two dimen¬ 
sions (Ghiglia and Pritt, 1998). At least two 


types of unwrapping algorithms have been ex¬ 
tended to three dimensions: branch cut algo¬ 


rithms (Salfity et ah, 2006; Huntley, 2001; Mark¬ 


lund et al. 

, 2007) and quality guided unwrapping 

(Cusack and Papadakis 

, 2002; 

Abdul- Rahman 


et ah, 2007, 2009). Branch cut algorithms iden¬ 


tify singularities in the phase volume and use 
them to identify and remove unwrapping paths 
that will lead to error. The phase is then 
unwrapped using a breadth-first spanning tree. 
Quality guided algorithms use some heuristic to 
evaluate the quality of each possible connection 
between adjacent points. The highest quality 


edges are unwrapped first. Examples of these 
quality heuristics are local variance of phase val¬ 
ues or of phase differences. This research in¬ 
volved a three dimensional phase volume, so I 
implemented these 3-D unwrapping algorithms 
and applied them to the Utah array data. 

We know from examining the voltage data 
that the LFP is highly correlated over the array 
channels. Since all the channels are oscillating 
together, we observe that the phase should differ 
by no more than 27r across the array at any given 
time. This observation gives us a heuristic for 
evaluating the performance of phase unwrapping 
algorithms: calculating the mean phase differ¬ 
ence across the array for the unwrapped phase. 
Though the performance of the algorithms varies 
somewhat with the passband chosen, the qual¬ 
ity guided algorithms generally outperformed the 
branch cut algorithm. For the rest of the anal¬ 
ysis, the data were unwrapped with the qual¬ 
ity guided unwrapping algorithm using the lo¬ 
cal variance of phase differences as the quality 
heuristic. 

5.3 Fitting the Model 

The planar wave model is linear. Therefore, 
modeling the unwrapped phase using the estab¬ 
lished technique of linear regression is appropri¬ 
ate. The model is the simple linear model 

4>{x,y,t) = l3o +^ix +^ 2 y + (10) 

However, we must be able to account for changes 
over time so we cannot fit the entire phase vol¬ 
ume simultaneously. Instead, we fit multiple 
separate models to overlapping time windows of 
data. To estimate the wave parameters at time 
t using a window size of Sms, we use data from 
time [t —2]ms to [t-h2]ms. This means that mod¬ 
eling a 2 second section of data requires fitting 
1,996 separate models, one for each time point 
from Sms to 1,998ms. 

Once all the models are fitted, we can esti¬ 
mate the wave parameters along with standard 
confidence intervals. We can also examine the 
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B? values as one indicator of whether the planar 
wave model is a good fit to the data. Finally, 
we can perform formal hypothesis tests to de¬ 
termine whether there is any statistically signif¬ 
icant linear dependence of phase on spatial posi¬ 
tion. If we fail to reject the null hypothesis that 
/3i = /32 = 0, we can conclude that we have no 
evidence for wavelike activity. 


By modeling narrow-band filtered data using 
3D phase unwrapping and linear regression we 
can produce interpretable estimates of wave pa¬ 
rameters on a single trial basis. Figure shows 
the results of this analysis for one trial. In this 
trial we see two primary wave propagation direc¬ 
tions and a switch between them. For most of 
the trial, the speed is around 40 cm/s, which is 


physiologically reasonable (Rubino et ah, 2006) 


At the time when the wavelike activity switches 
direction, the speed increasing drastically as the 
LFP activity becomes temporarily more syn¬ 
chronous. Additionally, during this transition we 
fail to reject the null hypothesis of the test de¬ 
scribed above. Overall, the null hypothesis was 
rejected for over 98% of the time points analyzed 
in the 17-18Hz band. This indicates that the 
wavelike activity is a highly persistent feature of 
the pre-motor neural system. 


5.4 Evaluating the Model 


The linear regression framework allows us to use 
standard analysis techniques to evaluate model 
performance. For instance, identifying outliers 
that result from isolated noise or unwrapping 
errors can explain some decreases in values 
such as that observed at around time 950 in fig¬ 
ure Figure shows the residuals for one of 
the regressions fit to data in this time period. 
The residuals reveal outliers limited to two chan¬ 
nels. Identifying and removing outlying channels 
could allow improved model fitting. 


CN 
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Fig. 10: The distribution over 6x and 6y 
used to generate artifieial data. 


6 Artificial Data 


In order to assess the improved and previous 
analysis techniques, we created artificial data 
that mimicked the bimodal wave direction dis¬ 
tribution observed in the real data. The ar¬ 
tificial data were generated using model pa¬ 
rameters generated using temporally smoothed 
Metropolis-Hastings sampling from a bimodal 
distribution. This distribution is shown in fig¬ 
ure The frequency term of the artificial data 


model is specified to be 17.5Hz. The wave pa¬ 
rameters of an example trial of artificial data is 
shown in figure [H By adding various types of 
noise to this artificial data, we can evaluate the 
performance of wave detection methodology. 

When no noise is added to the artificial data, 
both the original and improved methodologies 
do a good job of estimating both the speed and 
direction of the waves. Both methods make some 
significant error on less than 0.5% of the data 
during times when the speed and direction of 
the artificial waves are changing rapidly. 

When noise is added to the system, differences 
between the methods become apparent. Figure 
M shows the mean estimation errors for three 
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Speed! (cm/s) Direction (mdiane) R-squared 





Fig. 8: Results from analyzing a single trial of data in the 17-18Hz band. For 
eaeh plot, time runs along the horizontal axis and the values for eaeh time point 
are ealeulated using the linear regression eentered at that time point. The top plot 
shows the R? values. The middle plot shows estimated wave propagation direetion. 
The bottom plot shows estimated wave speed. The shaded areas are 95% eonfidenee 
intervals. Vertieal grey lines indieate the beginning and ending of the waiting period. 
Though the switeh in direetion in this trial oeeurs near the end of the waiting period, 
there is no elear pattern of this type apparent over all trials. 
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Fig. 9: Residuals plotted against model terms for a regression with low . The 
highlighted residuals are outlying data points from two ehannels. The sehematie in 
the eenter shows the loeation of these ehannels. 
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Fig. 11: An example of the direetion and 
speed values for an artifieial data trial. 



Speed 


different methods when filtered white noise is 
added to the artificial data. The three methods 
are (1) the original method from 
( 2006| ), (2) the original estimation 
ing the new narrow band filtering, and (3) the 
new method. The estimation errors from the 
three models show that the narrow band filtering 
is very important to both speed and direction es¬ 
timation while the new model fitting technique 
further improves estimation of speed. 
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7 Conclusions and Future 
Work 

The work reported in this paper makes a new 
scientific claim about the pattern of wavelike ac¬ 
tivity observed in monkey motor cortex: it is a 
persistent pattern, characteristic of the system. 
This contrasts with previous work which sug¬ 
gested that the wavelike activity was episodic. 
The methodology developed in this research al¬ 
lows single trial analysis of the data, and simu¬ 
lated data confirms that it results in improved 
estimation of wave parameters. This improve- 


Fig. 12: Mean error measures for estima¬ 
tion of artifieial data parameters. in- 

dieates the original method, ^^old-hfiltering’^ 
uses the original estimation teehnique with 
the new filtering teehnique, and ^^new’’ uses 
the entire improved proeessing pipeline. The 
apparently heavy tails on the box plots repre¬ 
sent less than half a pereent of the data but 
appear heavy beeause eaeh box plot is gener¬ 
ated using around 800,000 points. 
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ment in analysis methodology opens the door 
for additional investigation of the properties and 
purpose of traveling waves in the neural system. 

The largest estimation improvements resulted 
from modifications of the filtering procedure 
used before extracting instantaneous phase. Be¬ 
cause of the characteristics of the Hilbert trans¬ 
form and analytic signal, it is very important 
to apply a narrow band filter. A principled 
method for selecting the appropriate width and 
location of the pass band for filtering is a mat¬ 
ter for future investigation. It is possible to re¬ 
peat the analysis for multiple narrow frequency 
bands, thereby covering a wide frequency space. 
However, methodology for integrating the results 
from multiple frequency bands has not yet been 
developed. It is scientific question whether is is 
appropriate to treat LFP data as a truly narrow 
band signal. The current planar wave model has 
only a single term for frequency, so if we wish to 
model data with multiple frequency components 
the model would need to be extended. 


Phase-based modeling of wavelike activity is 
the most common approach, but it would also 
be possible to develop a method for modeling 
the LFP voltage directly. Since the voltage 
model for a traveling wave is nonlinear, fitting 
such a model would require nonlinear estima¬ 
tion techniques. However, modeling the volt¬ 
age would avoid the phase extraction and phase 
unwrapping process which introduces strict fil¬ 
tering requirements and possible unwrapping er¬ 
ror. If the planar wave model cannot be fit 
successfully to data, other voltage domain tech¬ 
niques are possible. Gabriel and Eckhorn (2003) 
present a method based on pairwise correla¬ 
tions between channels, successfully demonstrat¬ 
ing wavelike activity in monkey striate cortex us¬ 
ing a linear 15-electrode array. 


In summary, the work presented here con¬ 
tinues the task of creating an automated, sta¬ 
tistically principled detection methodology for 
the detection of planar waves in neural systems. 
This could aid the community of neuroscien¬ 
tists attempting to determine the prevalence and 


functional significance of such spatio-temporal 
patterns. 
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